import scipy.constants as C
import matplotlib

# matplotlib.use("PyQt6")
matplotlib.use("tkagg")
import matplotlib.pyplot as plt

plt.ion()
import scipy
import numpy
# import cst.results
# proj = cst.results.ProjectFile(r"C:\Users\18213\Documents\WeChat Files\wxid_awme81lgohso22\FileStorage\File\2024-06\test_SC.cst")


# ctrl insert
x = numpy.linspace(0, 10, 1000)
plt.figure()
plt.plot(x, scipy.special.j0(x), label="J0")
plt.plot(x, scipy.special.j1(x), label="J1")
plt.legend()

# plt.show(block=True)
N_zeros = 70#70
mu_0_ps = scipy.special.jn_zeros(0, N_zeros)


def Es(z, Q, a, b, ):
    arr = numpy.exp(-mu_0_ps * numpy.abs(z) / a) * (
            2 / mu_0_ps * scipy.special.j1(mu_0_ps * b / a) / scipy.special.j1(mu_0_ps)) ** 2 * numpy.sign(z)
    return Q / (2 * numpy.pi * C.epsilon_0 * b ** 2) * arr.sum()

vec_Es = numpy.vectorize(Es,)

z = numpy.linspace(-10e-3, 10*1e-3, 1000)
Q=1e-9
a=5*1e-3
b=4*1e-3

Es_ = vec_Es(z,Q,a,b)
import pandas
df = pandas.read_csv(r"C:\Users\18213\Desktop/2.csv",header= None)
plt.figure()
plt.plot(z*1e3,Es_ , label="Es, calculated")
plt.plot(df[0],df[1],'.-', label = "Es, CST")
plt.title("N= %d, Esmax = %.2e"%( N_zeros, Es_.max()))
plt.legend()
# print(Es_.max())